﻿using LxBasic;

using System;
using System.Diagnostics;
using System.Threading.Tasks;

namespace LxCore
{
    /*
 * 备注：
 * 在radAtDirs存储的rad值都是相对值，是减去了辐射背景强度的
 * 
 * 思路：
 * 1.FormBoundValue 一步只更新一次，计算出
 * 1.1 cell中的积分，按9.1，存到RadIntergral数组，对于无折射吸收的可不计算；
 * 1.2 边界中积分，按9.2，需要先由内部计算出边界网格上的值，积分已存储在数组1中，总体计算出等号右半部分后Copy到其余（漫反射应该是一样的）（因为边界温度可能已经改变-对流导致，所以必须在这里计算边界条件）
 * 2.FormAb step1和step2，每个方向计算一次
 * 3.UpdateBoundValue 一步只更新一次，需要更新出边界的能量源项，积分项存储到1，热源项存储到2/3/4
 * 
 * 存储：
 * RadAtDir[1] 边界存储9.2中积分项，即
 * RadAtDir[2] 边界存储热流a
 * RadAtDir[3]边界存储热流b
 * RadAtDir[4]边界存储热流
 */
    public class DOMRadiation : MatrixBase
    {
        public static readonly double StefanBoltzmann = 5.6703e-8;//黑体辐射常数sigma,W/(m2.K4)

        //protected double[] RadiaIntensInter;//存储的9.1 9.2的值，边界上除了积分，还包含壁面发射
        public BoundHead[] RadiationBounds;
        protected double BackgroundIntensity;//背景辐射强度
        public DiscreteRadiation[] DiscreteRadiations;
        public IField<double> AbsorbCoeff;//AbsorptionCoefficient 1/m;
        public IField<double> ScattCoeff;//ScatteringCoefficient 1/m;
        protected int CurrentDiscreteIndex;
        protected DomTypes DomRadType;
        //public DomToHeatSource HeatSourceOfRadiation;
        //public DomToHeatFlux HeatFluxOfRadiation;
        protected bool IsADLUExist = false;//ADLU不用每次都更新
        //share
        public ScalarBase Temperature;
        public void PreCheckCreate(MatrixBase arrayHolder)
        {
            Grad = arrayHolder.Grad;//不使用梯度，屏蔽
            InnX = arrayHolder.InnX;//防止被重复创建
            AInnO = arrayHolder.AInnO;
            AInnN = arrayHolder.AInnN;
            ADiag = arrayHolder.ADiag;
            B = arrayHolder.B;
        }
        public override void CheckCreate(FvZone zone)
        {
            base.CheckCreate(zone);
            //LxTool.CheckCreate(ref RadiaIntensInter, Cells.Length);//base.CheckCreate之前Cell未赋值
            LxBasic.CheckCreator.ArrayStrict(ref RadiationBounds, zone.Bounds.Count);
            //LxTool.CheckCreate(ref ScalarValue, out ScalarValueWithRef scalar);
            //scalar.SetScalar(zone, x);
            //LxTool.CheckCreate(ref HeatSourceOfRadiation);
            //LxTool.CheckCreate(ref HeatFluxOfRadiation);
        }

        public void SetZone_DomTerm(ScalarBase temperature, DomTypes domRadType, double backgroundTemperature)
        {//IFieldHost<double> absorbCoeff, IFieldHost<double> scattCoeff,
            AbsorbCoeff = ConstField.Double0;// absorbCoeff;
            ScattCoeff = ConstField.Double0; //scattCoeff;
            this.Temperature = temperature;
            BackgroundIntensity = StefanBoltzmann * Math.Pow(backgroundTemperature, 4) / Math.PI;
            //(ScalarValue as ScalarValueWithRef).SetRefValue(BackgroundIntensity);
            DomRadType = domRadType;
            int N = 0;
            switch (DomRadType)
            {
                case DomTypes.S4: N = 4; break;
                case DomTypes.S6: N = 6; break;
                case DomTypes.S8: N = 8; break;
                default: throw new NotImplementedException();
            }
            checkCreateDirs(N);
            setSnInOctant1(N);
            setOtherDirsByOctant1();

            //HeatSourceOfRadiation.SetParameter(Zone, DiscreteRadiations[1].x, DiscreteRadiations[2].x);
            //HeatFluxOfRadiation.SetParameter(Zone, DiscreteRadiations[3].x);
        }
        public void SetBound_FixRadiation(FvBound fvThread, IField<double> radationTemperature, IField<double> Emissivity,
            IField<double> Reflectivity, IField<double> Transmissivity, bool IsReflectivityAutoCal)//这里的radiation和边界条件一样，调用GetCellValue
        {
            LxBasic.CheckCreator.Instance(ref RadiationBounds[fvThread.ID], out FixRadiationBound b);
            b.SetRadiation(Emissivity, Reflectivity, Transmissivity, IsReflectivityAutoCal);
            b.SetRadiationTemperature(radationTemperature);
        }
        public void SetBound_FollowTemperature(FvBound fvThread, IField<double> Emissivity,
            IField<double> Reflectivity, IField<double> Transmissivity, bool IsReflectivityAutoCal)
        {
            LxBasic.CheckCreator.Instance(ref RadiationBounds[fvThread.ID], out FollowTempeartureBound b);
            b.SetRadiation(Emissivity, Reflectivity, Transmissivity, IsReflectivityAutoCal);
        }
        public void SetBound_Inactive(FvBound fvThread)
        {
            throw new NotImplementedException();
        }
        private void checkCreateDirs(int N)
        {
            throw new NotImplementedException();
            //int count = N * (N + 2);
            //if (DiscreteRadiations == null || DiscreteRadiations.Length != count)
            //{
            //    DiscreteRadiations = new DiscreteRadiation[count];
            //    for (int i1 = 0; i1 < count; i1++)
            //    {
            //        ref var disRad = ref DiscreteRadiations[i1];
            //        LxBasic.CheckCreator.ArrayStrict(ref disRad.x, Zone.Cells.Length);
            //        LxBasic.CheckCreator.ArrayStrict(ref disRad.ADiag, Zone.Cells.Length);
            //        LxBasic.CheckCreator.ArrayStrict(ref disRad.ALUO, Zone.Faces.Length);
            //        LxBasic.CheckCreator.ArrayStrict(ref disRad.ALUN, Zone.Faces.Length);
            //    }
            //    SetIsInitialized(false);
            //    //IsInitialized = false;
            //    IsADLUExist = false;
            //}
        }
        private void setSnInOctant1(int n)
        {
            int half_n = n / 2;
            double miu0;
            switch (n)
            {
                case 4:
                    miu0 = 0.2958759;
                    double s4w1 = Math.PI / 6;
                    DiscreteRadiations[0].Weight = s4w1;
                    DiscreteRadiations[1].Weight = s4w1;
                    DiscreteRadiations[2].Weight = s4w1;
                    break;
                case 6:
                    miu0 = 0.1838670;
                    double s6w1 = 0.1609517, s6w2 = 0.3626469;
                    DiscreteRadiations[0].Weight = s6w1;
                    DiscreteRadiations[1].Weight = s6w2;
                    DiscreteRadiations[2].Weight = s6w1;
                    DiscreteRadiations[3].Weight = s6w2;
                    DiscreteRadiations[4].Weight = s6w2;
                    DiscreteRadiations[5].Weight = s6w1;
                    break;
                case 8:
                    miu0 = 0.1422555;
                    double s8w1 = 0.1712359, s8w2 = 0.0992284, s8w3 = 0.4617179;
                    DiscreteRadiations[0].Weight = s8w1;
                    DiscreteRadiations[1].Weight = s8w2;
                    DiscreteRadiations[2].Weight = s8w2;
                    DiscreteRadiations[3].Weight = s8w1;
                    DiscreteRadiations[4].Weight = s8w2;
                    DiscreteRadiations[5].Weight = s8w3;
                    DiscreteRadiations[6].Weight = s8w2;
                    DiscreteRadiations[7].Weight = s8w2;
                    DiscreteRadiations[8].Weight = s8w2;
                    DiscreteRadiations[9].Weight = s8w1;
                    break;
                default:
                    throw new NotImplementedException();
            }
            double[] miu = new double[half_n];
            //计算cosine
            miu[0] = miu0;
            double miu0_pow2 = miu0 * miu0;
            for (int i = 1; i < half_n; i++)
            {
                miu[i] = Math.Sqrt(
                   miu0_pow2 + i * 2 * (1 - 3 * miu0_pow2) / (n - 2)
                    );
            }
            //计算第一象限的方向
            int radIndex = 0;
            for (int i = 0; i < half_n; i++)
            {
                int jcount = half_n - i;
                for (int j = 0; j < jcount; j++)
                {
                    DiscreteRadiations[radIndex].Direction = (new Vector3d(miu[i], miu[j], miu[half_n - 1 - i - j]));
                    radIndex++;
                }
            }
        }
        private void setOtherDirsByOctant1()//设置好第一象限后，更新其他方向
        {
            int groupCount = DiscreteRadiations.Length / 8;
            for (int i1 = 0; i1 < groupCount; i1++)
            {
                double weight = DiscreteRadiations[i1].Weight;
                Vector3d dir = DiscreteRadiations[i1].Direction;
                DiscreteRadiations[i1 + groupCount].SetDir(new Vector3d(-dir.X, dir.Y, dir.Z), weight);
                DiscreteRadiations[i1 + 2 * groupCount].SetDir(new Vector3d(dir.X, -dir.Y, dir.Z), weight);
                DiscreteRadiations[i1 + 3 * groupCount].SetDir(new Vector3d(dir.X, dir.Y, -dir.Z), weight);
            }
            groupCount = DiscreteRadiations.Length / 2;
            for (int i1 = 0; i1 < groupCount; i1++)
            {
                DiscreteRadiations[i1 + groupCount].SetDir(-DiscreteRadiations[i1].Direction, DiscreteRadiations[i1].Weight);
            }
        }

        public override void Initialize()
        {
            throw new NotImplementedException();
            //MultiThreads.RunSubThreads((threadInfo) =>
            //{
            //    (int i1, int endp1) = Zone.CellIdx.SliceToThread( threadInfo);
            //    for (int i2 = 0; i2 < DiscreteRadiations.Length; i2++)
            //    {
            //        Array.Clear(DiscreteRadiations[i2].x, i1, endp1 - i1);
            //    }
            //});
        }
        protected void formBoundValue(FvBound bound, IField<double> boundTemperature, IField<double> emissivity, IField<double> reflectivity, ThreadInfo threadInfo)
        {
            throw new NotImplementedException();
            //double[] radIntergral = DiscreteRadiations[0].x;
            //(int start, int endp1) = bound.FaceIdx.Split( threadInfo);//边界上是face统计的
            //var emi = emissivity.GetFaceVars(bound);
            //var boundT = boundTemperature.GetCellVars(bound);
            //var refl = reflectivity.GetFaceVars(bound);
            //for (int i1 = start; i1 < endp1; i1++)
            //{
            //    ref FvFace face = ref Faces[i1];
            //    double radInt = radIntergral[face.NSideCell];//取出上一步积分结果
            //    double boundradI = (emi.GetFaceVar(i1) * StefanBoltzmann * Math.Pow(boundT.GetCellVar(face.NSideCell), 4)
            //         + refl.GetFaceVar(i1) * radInt)
            //         / Math.PI - BackgroundIntensity;
            //    radIntergral[face.NSideCell] = boundradI;
            //}
            //int cellStart = Faces[start].NSideCell, length = endp1 - start;
            //for (int dirs = 1; dirs < DiscreteRadiations.Length; dirs++)//当前存储在0上，因此拷贝到其他
            //{
            //    Array.Copy(radIntergral, cellStart, DiscreteRadiations[dirs].x, cellStart, length);
            //}
        }
        //public override void FormBoundValue(ThreadInfo threadInfo)
        //{
        //    //for (int i1 = 0; i1 < RadiationBounds.Length; i1++)
        //    //{
        //    //    RadiationBounds[i1].FormBoundValue(this, Zone.Bounds[i1], threadInfo);
        //    //}
        //}
        public override void FormABOnFaces(ThreadInfo threadInfo)
        {
            throw new NotImplementedException();
            //Vector3d radDir = DiscreteRadiations[CurrentDiscreteIndex].Direction;
            //if (!IsADLUExist)
            //{
            //    (int i1, int endp1) = Zone.Inner.FaceIdx.Split(threadInfo);
            //    for (; i1 < endp1; i1++)
            //    {
            //        ref FvFace face = ref Faces[i1];
            //        double flux = Vector3d.DotProduct(face.Normal, radDir);
            //        if (flux >= 0)
            //        {
            //            ALUO[i1] = 0;
            //            ALUN[i1] = -flux;
            //        }
            //        else
            //        {
            //            ALUO[i1] = flux;
            //            ALUN[i1] = 0;
            //        }
            //    }
            //}
            //for (int i1 = 0; i1 < RadiationBounds.Length; i1++)
            //{
            //    RadiationBounds[i1].FormAbOnFaces(this, Zone.Bounds[i1], threadInfo);
            //}
        }
        protected void formAbOnFaces(FvBound bound, ThreadInfo threadInfo)
        {
            throw new NotImplementedException();
            //Vector3d radDir = DiscreteRadiations[CurrentDiscreteIndex].Direction;
            //(int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
            //for (; i1 < endp1; i1++)
            //{
            //    ref FvFace face = ref Faces[i1];
            //    double flux = Vector3d.DotProduct(face.Normal, radDir);
            //    if (!double.IsNegative(flux))
            //    {
            //        ALUO[i1] = 0;
            //        ALUN[i1] = -flux;
            //    }
            //    else
            //    {
            //        ALUO[i1] = flux * XInn[face.NSideCell];//x已经被赋予了边界条件//x[face.NSideCell] =
            //        ALUN[i1] = 0;
            //    }
            //}
        }
        public override void FormABOnCells(ThreadInfo threadInfo)
        {
            throw new NotImplementedException();
            //(int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
            //if (CurrentDiscreteIndex == 0)
            //{
            //    Array.Clear(B, i1, endp1 - i1);
            //}
            //if (IsADLUExist) return;
            //if (AbsorbCoeff.IsEqual0())
            //{
            //    for (; i1 < endp1; i1++)
            //    {
            //        ref FvCell cell = ref Cells[i1];
            //        ADiag[i1] = -cell.CalAllAdd(ALUN, ALUO);
            //    }
            //}
            //else
            //{
            //    var absorb = AbsorbCoeff.GetCellVars(Zone.Inner);
            //    for (; i1 < endp1; i1++)
            //    {
            //        ref FvCell cell = ref Cells[i1];
            //        ADiag[i1] = -cell.CalAllAdd(ALUN, ALUO) + absorb.GetCellVar(i1) * Cells[i1].Volumn;
            //    }
            //}
        }
        public override void UpdateBoundsValue(ThreadInfo threadInfo)
        {
            //if (!AbsorbCoeff.IsEqual0() || !ScattCoeff.IsEqual0())
            //{
            //    //这时候还需要对内部积分
            //    throw new NotImplementedException();
            //}
            //for (int i1 = 0; i1 < RadiationBounds.Length; i1++)
            //{
            //    RadiationBounds[i1].UpdateBoundValue(this, Zone.Bounds[i1], threadInfo);
            //}
            //if (threadInfo == 0 && !IsADLUExist)
            //{
            //    IsADLUExist = true;
            //}
        }
        protected void updateBoundValue(FvBound bound, IField<double> emissivity, ThreadInfo threadInfo)
        {
            throw new NotImplementedException();
            //double[] tempera = Temperature.XInn;
            //double[] radIntergral = DiscreteRadiations[0].x;
            //double[] qflux_a = DiscreteRadiations[1].x;
            //double[] qflux_b = DiscreteRadiations[2].x;
            //double[] qflux = DiscreteRadiations[3].x;
            //(int i1, int endp1) = bound.FaceIdx.Split(threadInfo);//边界上是face统计的
            //var emiss = emissivity.GetFaceVars(bound);
            //for (; i1 < endp1; i1++)
            //{
            //    ref FvFace face = ref Faces[i1];
            //    double radInt = radIntergral[face.NSideCell] = calBound_E_ByRadiationIntensityIntergral(ref face);//单位面积上的
            //    double emi = emiss.GetFaceVar(i1);
            //    double xishou = emi * radInt;//直接积分出来的是相对值
            //    double tem = tempera[face.NSideCell];
            //    double A = -emi * StefanBoltzmann;
            //    double flux = qflux[face.NSideCell] = A * Math.Pow(tem, 4) + xishou;
            //    qflux_a[face.NSideCell] = 4 * A * Math.Pow(tem, 3);
            //    qflux_b[face.NSideCell] = xishou - 3 * A * Math.Pow(tem, 4);
            //}
        }
        public void SetCurrentDir(int dirIndex)
        {
            if (dirIndex >= 0 && dirIndex < DiscreteRadiations.Length)
            {
                CurrentDiscreteIndex = dirIndex;
                ref var disRad = ref DiscreteRadiations[dirIndex];
                this.InnX = disRad.x;
                this.ADiag = disRad.ADiag;
                this.AInnO = disRad.ALUO;
                this.AInnN = disRad.ALUN;
            }
            else
            {
                throw new Exception("DOM辐射计算，方向设置出错");
            }
        }
        protected double calRadiationIntegral_inner(int cellID)
        {
            double sum = 0;
            foreach (var r in DiscreteRadiations)
            {
                sum += r.x[cellID] * r.Weight;
            }
            return sum / (4 * Math.PI);
        }
        protected double calBound_E_ByRadiationIntensityIntergral(ref FvInnFace face)
        {
            double sum = 0;
            int half = DiscreteRadiations.Length / 2;
            double scale = 0;
            for (int i1 = 0; i1 < half; i1++)
            {
                ref var rad1 = ref DiscreteRadiations[i1];
                double temp = Vector3d.DotProduct(face.Normal, rad1.Direction);
                if (!double.IsNegative(temp))
                {
                    temp *= rad1.Weight;
                    scale += temp;
                    sum += (rad1.x[face.OSideCell] + BackgroundIntensity) * temp;//还需要除以面积，后面统一除
                }
                else
                {
                    ref var rad2 = ref DiscreteRadiations[i1 + half];
                    temp *= -rad2.Weight;
                    scale += temp;
                    sum += (rad2.x[face.OSideCell] + BackgroundIntensity) * temp;//还需要除以面积，后面统一除
                }
            }
            return sum / face.Area;
        }
        public class FollowTempeartureBound : BoundHead
        {
            public IField<double> Emissivity;//发射率epsilon，=吸收率alpha
            public IField<double> Reflectivity;//反射率rho，建议设置自动计算
            public IField<double> Transmissivity;//透射率gamma
            public void SetRadiation(IField<double> Emissivity, IField<double> Reflectivity, IField<double> Transmissivity, bool IsReflectivityAutoCal)
            {
                this.Emissivity = Emissivity;
                this.Transmissivity = Transmissivity;
                if (IsReflectivityAutoCal)
                {
                    if (Emissivity is ParameterBase<double> emi && Transmissivity is ParameterBase<double> trans)
                    {
                        this.Reflectivity = new ConstField.Double(1 - emi.GetPara() - trans.GetPara());
                    }
                    else
                    {
                        throw new NotImplementedException();
                    }
                }
                else
                {
                    this.Reflectivity = Reflectivity;
                }
            }
            //public override void FormBoundValue(object inMatrix, FvBound fvThread, ThreadInfo threadInfo)
            //{
            //    throw new NotImplementedException();
            //    //var m = inMatrix as DOMRadiation;
            //    //m.formBoundValue(fvThread,new CellArrayScalar() { CommonArray = m.Temperature.XInn }, Emissivity, Reflectivity, threadInfo);
            //}
            //public override void FormAbOnFaces(object inMatrix, FvBound fvThread, ThreadInfo threadInfo)
            //{
            //    var m = inMatrix as DOMRadiation;
            //    m.formAbOnFaces(fvThread, threadInfo);
            //}
            //public override void UpdateBoundValue(object inMatrix, FvBound fvThread, ThreadInfo threadInfo)
            //{
            //    var m = inMatrix as DOMRadiation;
            //    m.updateBoundValue(fvThread, Emissivity, threadInfo);
            //}
        }
        public class FixRadiationBound : FollowTempeartureBound
        {
            public IField<double> RadiationTemperature;
            public void SetRadiationTemperature(IField<double> radiationTemperature)
            {
                this.RadiationTemperature = radiationTemperature;
            }
            //public override void FormBoundValue(object inMatrix, FvBound fvThread, ThreadInfo threadInfo)
            //{
            //    var m = inMatrix as DOMRadiation;
            //    m.formBoundValue(fvThread, RadiationTemperature, Emissivity, Reflectivity, threadInfo);
            //}
        }
        public class InactiveBound
        {

        }
        public struct DiscreteRadiation
        {
            public double[] ADiag;
            public double[] ALUO;
            public double[] ALUN;
            public double[] x;
            public double Weight;
            public Vector3d Direction;
            public void SetDir(Vector3d dir, double weight)
            {
                this.Direction = dir;
                this.Weight = weight;
            }
        }

    }
    public enum DomTypes : byte
    {
        S4, S6, S8
    }
}
//public class DomToHeatSource : Value2Base<double>
//{
//    double[] aArray;
//    double[] bArray;
//    int deltaID_cell_face;
//    public void SetParameter(FvZone zone, double[] a, double[] b)
//    {
//        deltaID_cell_face = zone.InnerIndex.CellIDs.EndP1 - zone.InnerIndex.FaceIDs.EndP1;
//        aArray = a;
//        bArray = b;
//    }
//    public override void GetFaceValue(out double S_Ad, out double S_b, int faceIndex)
//    {
//        int cellid = faceIndex + deltaID_cell_face;
//        S_Ad = aArray[cellid];
//        S_b = bArray[cellid];
//    }
//}
//public class DomToHeatFlux : IFieldHost<double>
//{
//    double[] sArray;
//    int deltaID_cell_face;
//    public void SetParameter(FvZone zone, double[] s)
//    {
//        deltaID_cell_face = zone.InnerIndex.CellIDs.EndP1 - zone.InnerIndex.FaceIDs.EndP1;
//        sArray = s;
//    }
//    public override double GetFaceValue(int faceIndex)
//    {
//        return sArray[faceIndex + deltaID_cell_face];
//    }
//}
//public void SetBound_FixTemperature(FvBound fvThread, IFieldHost<double> Emissivity,
//    IFieldHost<double> Reflectivity, IFieldHost<double> Transmissivity, bool IsReflectivityAutoCal)
//{
//    //LxTool.CheckCreate(ref RadiationBounds[fvThread.ID], out RadBound b);
//    //b.SetHead(fvThread);
//    //b.SetRadiation(Emissivity, Reflectivity, Transmissivity, IsReflectivityAutoCal);
//}
//double[] radIntergral = DiscreteRadiations[0].x;
//double[] qflux_a = DiscreteRadiations[1].x;
//double[] qflux_b = DiscreteRadiations[2].x;
//double[] qflux = DiscreteRadiations[3].x;
//foreach (var b in RadiationBounds)
//{
//    var emiss = b.Emissivity;
//    b.BoundIndex.FaceIDs.GetSeg(  threadInfo);//边界上是face统计的
//    for (; i1 < endp1; i1++)
//    {
//        ref FvFace face = ref Faces[i1];
//        double radInt = radIntergral[face.NSideCell] = calBound_E_ByRadiationIntensityIntergral(ref face);//单位面积上的
//        double emi = emiss.GetFaceValue(i1);
//        double xishou = emi * radInt;//直接积分出来的是相对值
//        double tem = tempera[face.NSideCell];
//        double A = -emi * StefanBoltzmann;
//        double flux = qflux[face.NSideCell] = A * Math.Pow(tem, 4) + xishou;
//        qflux_a[face.NSideCell] = 4 * A * Math.Pow(tem, 3);
//        qflux_b[face.NSideCell] = xishou - 3 * A * Math.Pow(tem, 4);
//    }
//}
//foreach (var b in RadiationBounds)
//{
//    b.BoundIndex.FaceIDs.GetSeg(  threadInfo);
//    for (; i1 < endp1; i1++)
//    {
//        ref FvFace face = ref Faces[i1];
//        double flux = Vector3d.DotProduct(face.Normal, radDir);
//        if (!double.IsNegative(flux))
//        {
//            ALUO[i1] = 0;
//            ALUN[i1] = -flux;
//        }
//        else
//        {
//            ALUO[i1] = flux * x[face.NSideCell];//x已经被赋予了边界条件//x[face.NSideCell] =
//            ALUN[i1] = 0;
//        }
//    }
//}
//public ScalarValueWithRef RadiationValue;
//protected RadBound[] RadiationBounds;
//public class RadBound : BoundBase
//{
//    public void SetRadiation(IFieldHost<double> Emissivity, IFieldHost<double> Reflectivity, IFieldHost<double> Transmissivity, bool IsReflectivityAutoCal)
//    {
//        this.Emissivity = Emissivity;
//        this.Transmissivity = Transmissivity;
//        if (IsReflectivityAutoCal)
//        {
//            if (Emissivity is ConstField.Double emi && Transmissivity is ConstField.Double trans)
//            {
//                this.Reflectivity = new ConstField.Double(1 - emi.Number - trans.Number);
//            }
//            else
//            {
//                throw new NotImplementedException();
//            }
//        }
//        else
//        {
//            this.Reflectivity = Reflectivity;
//        }
//    }
//    public IFieldHost<double> Emissivity;//发射率epsilon，=吸收率alpha
//    public IFieldHost<double> Reflectivity;//反射率rho，建议设置自动计算
//    public IFieldHost<double> Transmissivity;//透射率gamma
//}

//double[] tempera = Temperature.x;
//double[] radIntergral = DiscreteRadiations[0].x;
//foreach (var b in RadiationBounds)
//{
//    var emiss = b.Emissivity; var reflect = b.Reflectivity;
//    b.BoundIndex.FaceIDs.GetSeg(out int start, out int endp1, threadInfo);//边界上是face统计的
//    for (int i1 = start; i1 < endp1; i1++)
//    {
//        ref FvFace face = ref Faces[i1];
//        double radInt = radIntergral[face.NSideCell];//取出上一步积分结果
//        double boundradI = (emiss.GetFaceValue(i1) * StefanBoltzmann * Math.Pow(tempera[face.NSideCell], 4)
//             + reflect.GetFaceValue(i1) * radInt)
//             / Math.PI - BackgroundIntensity;
//        radIntergral[face.NSideCell] = boundradI;
//    }
//    int cellStart = Faces[start].NSideCell, length = endp1 - start;
//    for (int dirs = 1; dirs < DiscreteRadiations.Length; dirs++)//当前存储在0上，因此拷贝到其他
//    {
//        Array.Copy(radIntergral, cellStart, DiscreteRadiations[dirs].x, cellStart, length);
//    }
//}
//switch (DomRadType)
//{
//    case DomTypes.S4:
//        checkCreateDirs(4 * 6);
//        double x1 = 0.2958759, x2 = 0.9082483, w1 = Math.PI / 6;
//        RadiationAtDirs[0].SetDir(new Vector3d(x1, x1, x2), w1);
//        RadiationAtDirs[1].SetDir(new Vector3d(x1, x2, x1), w1);
//        RadiationAtDirs[2].SetDir(new Vector3d(x2, x1, x1), w1);
//        break;
//    case DomTypes.S6:
//        checkCreateDirs(6 * 8);
//        break;
//    case DomTypes.S8:
//        checkCreateDirs(8 * 10);
//        break;
//    default:
//        throw new NotImplementedException();
//}
///内部和边界的定义是不一样的，要注意
///内部就是公式中辐射强度求和
///边界处是边界条件
//public void CalRadiationIntensityIntergral(ThreadInfo threadInfo)
//{
//    double[] tempera = Temperature.x;
//    if (!AbsorbCoeff.IsEqual0() || !ScattCoeff.IsEqual0())
//    {
//        //这时候还需要对内部积分
//        throw new NotImplementedException();
//    }
//    foreach (var b in RadBounds)
//    {
//        var emiss = b.Emissivity; var reflect = b.Reflectivity;
//        b.ThreadIndex.FaceIDs.GetSeg(  threadInfo);//边界上是face统计的
//        for (; i1 < endp1; i1++)
//        {
//            RadiaIntensInter[i1] = calRadiationIntergral_bound(i1);
//            //RadiaIntensInter[i1] = (emiss.GetFaceValue(i1) * StefanBoltzmann * Math.Pow(tempera[i1], 4) - ReferenceRadiationIntensity)
//            //    + reflect.GetFaceValue(i1) * calRadiationIntergral_bound(i1);
//        }
//    }
//}
//public void CalBoundRadiationCondition(ThreadInfo threadInfo)
//{
//    double[] tempera = Temperature.x;
//    foreach (var b in RadBounds)
//    {
//        var emiss = b.Emissivity; var reflect = b.Reflectivity;
//        b.ThreadIndex.FaceIDs.GetSeg(  threadInfo);//边界上是face统计的
//        for (; i1 < endp1; i1++)
//        {
//            RadiaIntensInter[i1] = emiss.GetFaceValue(i1) * StefanBoltzmann * Math.Pow(tempera[i1], 4)
//                - ReferenceRadiationIntensity
//                + reflect.GetFaceValue(i1) * RadiaIntensInter[i1];
//        }
//    }
//}